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Abstract 


A two-dimensional, two-temperature, single fluid magnetohydrodynamics code which incorporates classical 
plasma transport coefficients and Hall effects has been developed to predict steady-state, self-field MPD 
thruster performance. The governing equations and numerical methods of solution are outlined and discussed. 
Experimental comparisons are used to validate model predictions. The model accurately predicts thrust and 
reproduces trends in the discharge voltage for discharge currents below experimentally measured onset values. 
However, because the model does not include electrode effects the calculated voltage drops are significantly 
lower than experimentally measured values. Predictions of thrust and flow efficiency are made for a matrix 
of fifteen cylindrical thruster geometries assuming a fully ionized argon propellant. A maximum predicted 
specific impulse of 1680 s is obtained for a thruster with an anode radius of 2.5 cm, a cathode radius of 0.5 
cm, and equal electrode lengths of 2.5 cm. A scaling relation is developed to predict, within limits, the onset 
of cylindrical, self-field thruster instability as a function of geometry and operating condition. 


Nomenclature 


A 

plasma ionization function 

^0 ,(ea:it) 

inlet (exit) velocity, m/s 

A c 

thruster exit area, m 2 

V 

discharge voltage, V 

B 

magnetic field, T 

v^t 

external volume, m 3 

e 

electric charge, C 

Veot 

total volume, m 3 

E 

electric field, V/m 

On 

finite difference coeffs. 

F 

thrust, N 

Tn 

differential equation coeffs. 

g 

gravitational acceleration, 9.8 m/s 2 

pn 
A m 

differential equation coeffs. 

hp 

specific impulse, s 

fo 

permittivity of free space, 8.854 x 10“ 12 C 2 /N-: 

j 

current density, A/m 2 

1 

viscosity, kg/m-s 

J 

discharge current, A 

Vf 

plasma flow efficiency 

la.c 

length (anode, cathode), m 

«t(e) 

ion (electron) thermal conductivity, Watt/m-K 

k D 

Boltzmann const., J/K 

A 

second coeff. of viscosity, kg/m-s 

m, M 

electron, ion mass, kg 

In A*( e ) 

ion (electron) Coulomb logarithm 

m 

mass flow rate, kg/s 

Mo 

permeability of free space 

n 

number density, m“ 3 

Vex 

electron-ion collision frequency, s'" 1 

P(i, c ) 

pressure (ion, electron), Pa 

( 

density function 

P exit 

exit plane pressure, Pa 

p 

mass density, kg/m 3 

Ptanfc 

background gas pressure, Pa 

<J 

electrical conductivity, mho/m 

P 

power, W 

r e 

electron collision time, s 

AQci 

energy exchange term, J/s 


electric potential, V 

r (a,c) 

radius (anode, cathode), m 


viscous dissipation function 

R 

gas constant, J/K-mol 


magnetic field streamline 

s 

general source term 

* 

viscous force vector 

S tC 

thrust chamber inner surface, m 2 

w 

convergence factor 


ion (electron) temperature, K 

Vce 

electron cyclotron frequency, s“ 1 

V 

velocity, m/s 

n 

electron Hall parameter 


Subscripts: 

i,e 

ion, electron 



r.0,2 

radial coordinates 
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I, Introduction 


The magnetoplasmadynamic (MPD) thruster has been advocated for a variety of propulsion applications, 
from near-Earth orbit raising maneuvers 1 to the long duration interplanetary missions envisioned by the 
National Space Exploration Initiative. 2 In its basic form, the MPD thruster consists of a cylindrical cathode 
surrounded by a concentric anode (Figure 1). An arc struck between the electrodes ionizes a gaseous 
propellant, and the interaction of the current with the self-induced magnetic field accelerates the plasma to 
produce thrust. Steady-state MPD thrusters have been operated at power levels approaching 600 kW, while 
pulsed, quasi-steady devices have operated in the megawatt power range. 3 The engine is robust and designed 
to provide low, continuous thrust at specific impulse values between 1,000 and 10,000 s. 

MPD thruster performance is currently limited by low thrust efficiency in the operating regimes of interest, 3,4,5 
with significant fractions of the applied power deposited into the anode. For low power steady-state thrusters, 
the anode power fraction may reach 80%, leaving only 20% for plasma ionization and acceleration. 3 Experi- 
mentally, the use of applied magnetic fields 3,4 and flared electrode geometries 6 have been shown to increase 
the specific impulse and improve thruster efficiency, although for reasons poorly understood at present. 
The specific impulse for a self-field MPD thruster is related to J 2 /rh, a parameter which is often used to 
characterize MPD thruster performance. 3 High values of J 2 /m correspond to predominantly electromagnetic 
acceleration, and provide higher values of specific impulse; low values of J 2 jrn correspond to predominantly 
electrothermal acceleration, and lower values of specific impulse. MPD thruster efficiency typically increases 
with increasing J 2 /m, but is limited by the occurrence of discharge voltage oscillations which result in severe 
electrode erosion. Several mechanisms have been proposed to explain the onset of these instabilities, includ- 
ing anode mass starvation, 7 flow choking due to enhanced back-EMF, 8 * 9 and the triggering of electrothermal 
and gradient-driven instabilities as the plasma approaches full ionization. 10,11,12 It appears that different 
operating conditions may trigger one or a combination of the proposed mechanisms, with a concomitant 
reduction in thruster performance and lifetime. 

MPD thrusters have undergone extensive experimental development since the 1960’s, 3 ' 5 but a comprehensive 
theoretical analysis has been hampered by the complex nature of the coupled electromagnetic and gasdynamic 
acceleration processes. Analytic and numerical simulations using ID and quasi-ID approximations of the 
magnetohydrodynamic fluid equations have provided valuable insights into both self-field and applied-field 
MPD thruster operation, but they are by definition constrained in their ability to predict detailed thruster 
performance. The removal of such limitations via 2D and quasi-2D modeling has become practical with the 
emergence of high-speed computational facilities, allowing model validation and refinement using the existing 
experimental data base, while in return establishing a theoretical basis to guide further experimentation. A 
review of MPD thruster technology and a survey of recent MPD thruster models is provided by Myers ei 
ai 3 

Self-Field MPD Thruster Scaling Studies. 


Steady improvements in MPD thruster performance have occurred primarily via the relatively expensive 
method of device fabrication and testing, with the more robust designs surviving this evolutionary process 
to undergo further testing and refinement. However, with the variety of thruster geometries currently 
under development, progress toward understanding the impact of geometric scaling on performance has been 
delayed by an incomplete empirical data base and by the replication of experimental effort. The limited 
research efforts which address performance variations due to geometric changes in the electrode shapes are 
reviewed below. 

One of the first comprehensive attempts to verify the effects of geometric variation on MPD thruster perfor- 
mance was reported by King, 13 who combined an experimental approach with a simplified ID gasdynamic 
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model to evaluate anode modifications in the Princeton benchmark thruster. King found that increased 
electrode lengths and decreased anode radii allow an increase in the discharge current before the onset of 
thruster instability, with the decreased anode orifice contributing most strongly to the increased stability. 
The results led to the design of a modified flared anode thruster, which consisted of two straight anode 
channels connected by a flared anode region. The design proved to have equivalent voltage-current charac- 
teristics and a slightly higher thrust efficiency compared to a straight cylindrical MPD thruster of similar 
dimensions. 

In a detailed experimental study, Gilland 14 compared the performance of the Princeton standard benchmark 
thruster and the modified flared anode thruster with their geometric half-scale counterparts. In each case, 
thruster performance was determined to be primarily a function of propellant mass flow rate and the ratio of 
the electrode radii, but was generally independent of the geometric scale of the thruster. The scaled thruster 
performance was similar for similar mass flow rates, although the half scale thrusters were less efficient at 
equal power levels and exhaust velocities. 

Schmidt 15 combined an empirical approach with a simple theoretical model to analyze the full scale and 
half scale Princeton benchmark and flared anode thrusters. Empirical data were employed to generate an 
analytic expression for voltage scaling in the devices. An expression was assumed for the current density 
distribution and used to evaluate thruster force components. An analytic expression was used to estimate 
the contributions by electrothermal forces, and the model was employed to analyze the performance of the 
full scale and half scale thrusters. The half scale thrusters, operated at half the power levels of the full scale 
versions, were found to have similar efficiencies and increased specific impulse, in general agreement with 
the experimental findings of Gilland. The model was to be used in conjunction with an Air Force Rocket 
Propulsion Laboratory (now Air Force Phillips Laboratory) experimental program to develop a variable 
geometry MPD thruster, designed specifically to address geometric scaling issues. Cassady et al 1G provide 
a brief description of the program, but little information has been published in the intervening years. 

Heimerdinger 17 performed an experimental evaluation of two thruster geometries, one a straight cylindrical 
channel and the second a straight cylindrical channel with diverging area near the exit plane. The constant 
area channel suffered from strong current concentrations near the inlet and exit regions, resulting in enhanced 
electrode erosion. The diverging electrode geometry diminished the discharge current densities and reduced 
the electrode erosion, and demonstrated that electrode geometries could be tailored to minimize ohmic 
heating and associated electrode erosion. 

Martinez-Sanchez 18 developed an analytic model using steady-state, quasi-ID MHD flow equations, and eval- 
uated a variety of electrode geometries. Optimum performance was obtained with a converging-diverging 
anode, which allowed more uniform current concentration along the electrodes. In descending order of per- 
formance were geometries which employed a diverging anode, a converging anode, and a straight cylindrical 
anode. The poor performance of the cylindrical electrode geometry was due to thermal pressure effects 
retarding the flow near the exit plane of the thruster. The results confirmed that electrode contouring can 
be effective in controlling the current distribution and associated energy dissipation effects. 

The experimental results of Kunii et al. 1 ^ confirmed that straight cylindrical thruster geometries did in- 
deed suffer from current concentrations near the base and tip of the cathode, with subsequent severe cath- 
ode erosion. A diverging anode was found to produce the most uniform current distribution; however, a 
converging-diverging anode showed severe current c onc entrations and electrode erosion, in contrast to the 
analytic models discussed above. In a comparison of cylindrical and flared thruster geometries, Uematsu 
et al . 20 demonstrated that the performance of a straight cylindrical thruster was improved by recessing the 
cathode, at the expense of increasing the cathode erosion rate. The flared anode thrusters were again shown 
to provide better performance, with better stability and less erosion, than cylindrical electrode geometries. 
In a recent modeling effort, Lefever-Button and Subramaniam 21 employed a quasi- Id model with finite rate 
kinetics to compare cylindrical and flared thruster geometries. For fixed operating conditions, the flared ge- 
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ometries produced higher exhaust velocities, reduced the current density concentrations near the exit plane, 
and increased the propellant ionization fraction compared to straight cylindrical thruster geometries. 

Preble 10 compiled a comprehensive data base of experimentally determined onset values for self-field, cylin- 
drical MPD thrusters operated with argon propellants. A stability model, developed to simulate the onset 
of electrothermal instabilities, was combined with a one-dimensional MHD code and used to evaluate onset 
conditions for a variety of thruster scale lengths. Through an extensive literature search, Preble catalogued 
several mechanisms which could lead to improved thruster performance. Longer electrode lengths, smaller 
interelectrode spacing, and propellant injection near the anode were found to significantly increase the value 
of J 2 /rn before the onset of thruster instability, allowing thrusters to operate more efficiently with higher spe- 
cific impulse. The trends predicted by the model correlated very well with experimentally determined onset 
conditions, and demonstrated that an electrothermal instability could precipitate global thruster instability 
over the range of geometric scales considered in the study. 

In addition to the onset of global thruster instability, plasma microinstabilities may arise which lead to an 
increase in the plasma resistivity, raising the voltage drop across the plasma and lowering thruster efficiency. 
Niewood et al. n combined a quasi- ID flow model with nonlinear dispersion relations to show that modified 
two-stream and electrothermal instabilities may arise in MPD thruster plasmas, with the onset of operational 
instability closely coupled to the occurrence of the electrothermal instability. Choueiri 22 and Choueiri et 
a/. 23 derived general nonlinear dispersion relations for current-driven plasma instabilities, and determined 
that lower-hybrid waves are the dominant microinstability for most MPD thruster operating regimes. The 
existence of lower-hybrid waves has been experimentally verified in both low power and high power MPD 
thrusters. 22,24 

Although limited in extent, significant information has been obtained from the experimental and numerical 
scaling studies outlined above. In general, flared electrode geometries provide more uniform current distri- 
butions, reduce electrode erosion, and enhance thruster stability. Straight cylindrical thrusters appear to 
operate more stably with longer electrodes, smaller interelectrode separations, and propellant injection near 
the anode. Half-scale reductions in thruster size might not significantly alter thruster performance for equal 
mass flow rates, although half-scale geometries are less efficient at equal power levels and exhaust velocities. 
Unfortunately, this impressive list of maxims is not yet sufficient to guide thruster designs toward an optimal 
configuration. Several uncertainties remain, such as the optimum scaling of electrode length versus thruster 
diameter, the impact of electrode length on thruster performance and efficiency, and the combined effect 
of geometric scale and operating condition on MPD thruster stability. The two-temperature MPD thruster 
model described below has been employed in an effort to address some of these lingering issues for cylindrical, 
self-field MPD thrusters. 

II. Two-Temperature MPD Thruster Model 

The two-temperature MPD thruster model is an extension of the single fluid, single temperature MPD 
code discussed in Reference [25]. Separate electron and ion energy equations have been incorporated into 
the steady-state, viscous, compressible, single fluid magnetohydrodynamic equations, which are written in 
cylindrical coordinates with assumed symmetry about the centerline. The plasma is assumed to be fully, 
singly ionized, and is described by a perfect gas equation of state. The inclusion of separate electron and ion 
temperatures permits refined calculations of the classical plasma transport coefficients, which provide more 
accurate estimates of the plasma voltage drop, viscous flow losses, and thruster efficiency. 

Ila. Electromagnetic Field Equations. 

The basic set of electromagnetic equations includes the full complement of Maxwell’s equations: 

(o) V • B = 0 (i) Vxfl=^j (1) 
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dB 


(c) V • E = p e /e o « 0 

which incorporate the steady-state, single fluid plasma approximations. A general Ohm’s law of the form: 


w Vx *=-7>r = 0 


j = v[E+(v*B)] 


( 2 ) 


is used to relate the current density j to the plasma velocity v and the electric ( E ) and magnetic ( B ) fields. 
Ion slip terms are neglected due to the assumption of full ionization. The electrical conductivity a is given 
by the classical Spitzer-Harm conductivity for a fully ionized plasma: 20 

j-3/2 


cr = 1.53 x 10" 


In A- 


where T e is the electron temperature in degrees- Kelvin, lnA e is the electron Coulomb logarithm: 

1.22 x 10 3 n 1 / 2 ] 


In A e « 23 - In 


n 3/2 


( 3 ) 


( 4 ) 


and n is the (single fluid) plasma number density expressed in particles per cubic meter. The electron Hall 
parameter $1 is the product of the electron cyclotron frequency (u; ce ) and the electron collision time (t c ): 


ft — W cc T e 


■ C) ( 


ma \ 
ne*J 


( 5 ) 


= 9.6 x 10 


16 


( T^B_\ 

^ZnlnA J 


where m is the electron mass, e is the electron charge, and B is the magnitude of the local magnetic field. 


Equation 2 may be rearranged to solve for the electric field distribution: 

1 


E = ~V<p = 


whose components are given by: 

E r = ~ 


d A 

dr 


J+~(jxB)j-(vxB) 

Ji + -g (jeB z — ji /lojj — (v e B z — v z Bg) 


( 6 ) 


F - d * ~ 0 
Ee ~ ~d0 = ° 


( 7 ) 


E z = - 


dz 


n 

jz + -g {j r Be - je B r ) 


- (v r B e - v 0 B T ) 


The radial electric field is integrated from cathode to anode to find the potential drop <f> across the plasma. 
The divergence of Eq. 6 may be used to generate a Poisson equation for the potential, or assuming quasineu- 
trality a simple Laplace equation may be solved to give the potential distribution, with the calculated plasma 
fall providing the necessary boundary conditions. The first method is rigorously correct, while the second 
method is computationally faster; both have been tested, and yield essentially identical results for the plasma 
potential distribution. 

A magnetic transport equation may be derived by combining Maxwell’s equations with the generalized Ohm’s 
law: 


V x E = - 


dB 

dt 


= V x 


V x B 
<r 


+ V x 


n 


lHocrB 


((V x B) x i?) j - V x (if x i?) = 0 


( 8 ) 
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(9) 


Defining (3 = Q/(aB) and ip = rB#, the azimuthal component of the above equation can be written: 


d 2 ip dip cty 3 2 t/> _ 

3r 2 + 71 dr +l2 dz + 73 ^ + dz 2 


where 


,'l 13a 

7i = - - + —~ET 

\r cr or 


72 


73 = 


, cr 3/3 \ 

~dz + ^ oC ™ r J 
1 3a 2/3a a 3/3 

■“ « a” + ^““5- + Mo<™ 

(j 02 r 3r 


0 


( 10 ) 


Ml) <7 


dv r v T dv z 

3T-7 + 37 


]) 


where the source term S incorporates effects of an applied magnetic field, and is equal to zero for the self- 
field thruster under consideration. Once the magnetic transport equation is solved for the value of Bq is 
readily obtained. The azimuthal magnetic field is a function of the total discharge current J, and through 
the relation ip — tBq ~ r(J/r) ~ J, the function ip(r t z) — constant may be used to represent lines of total 
enclosed current. 


Returning to the Maxwell equations, the radial and axial current densities are obtained from the calculated 
magnetic field distributions: 


Jr 


Jz 


1 dBe 
p o dz 



(ii) 


The azimuthal current density jo vanishes due to the absence of applied radial or axial magnetic fields in 
the self-field MPD thruster. 


lib. Fluid Equations. 

The fluid equations are based on the conservation equations of mass, momentum, and energy. Conservation 
of mass for a compressible fluid is given by: 

V • (pv) = 0 (12) 


which can be written: 


d(rpv T ) d(rpv x ) _ 
dr + dz 


(13) 


where p is the plasma mass density. Defining ( — rp, the mass conservation equation takes the simplified 
form: 


d(v r divz 
3r dz 


(14) 


which is solved for from which the value of p is determined. 

Momen tum. The conservation of momentum, including viscosity, is expressed in vector form as: 

/?(i»-V)f=-Vp+(jxij-|-$ (15) 
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(16) 


where p is the plasma pressure and 'P is the viscous force vector: 

= -|v [r?(v • V)] + r, [V 2 v + V (V ■ 6 )] + 2(Vr,-V)v + Vtj x (V x iT) 

o 


The viscosity 77 for a fully ionized plasma is given by: 


.27 


1.25 x 10 


_ 19 5 ( MmkpT^ 12 ^2 koT^J 


Sn'A V 


(17) 


where 7} is the ion temperature, M is the ion mass, m is the electron mass, e is the electron charge, k D is 
Boltzmann's constant, n‘=:l for a singly ionized plasma, and, for a fully ionized gas (a = 1), 


A = 2 


ln ( 1 + « 2 )“ 


= 0.386 


The radial component of the momentum equation is used to solve for the radial velocity v r : 

/ a n 2 drj\ 4 dv T / 77 dp\ 

' ) 3 dr \ r dr ) 


( dVr 1 


dv r 

!h 


+ -jzB„ 

or 


(A r, 

~ Vr \3r* 


2 ^ 3 r dr d 


^ dv r dr] irjd 2 v r d 2 v r 77 d 2 v z dp dv t 2 dp dv z 

+ dz dz f 3 <9r 2 ~ + 71 dz 2 + 3 drdz + dz ~dr ~ 3 dr dz 


and may be recast in the form: 

d 2 v r 


where 


and 


dr 2 


+ r 11 -^ + r 2 — 

' 1 r n_ T 1 r 


il 9Vjr 

dr 


3 3r, d 2 v r _ 

0z _1>r + ~4~W-~ St 


rj 

r 2 

r? 


( V dn _ 3/nv ^ 

\r dr A ) 


3 {dr, _ \ 

4 PVl ) 

(jl + IOr) 

\r 2 2 rdrj 


^ r 4 f dr 3 drdz dz dr 


r, d 2 v z dr, dv z 2 dr, dv z pv% 


3 dr dz 


+ 


(18) 


(19) 


(20) 


( 21 ) 


There are no radial or axial magnetic fields in the self-field MPD thruster, hence no azimuthal momentum 
is generated. The axial momentum component is given by: 


( dv z dv z \ 

f T'ar + ’’-a7j 


+ 


3 dz 2 

which yields an equation for the axial velocity v z \ 


f p dp\ dv 2 4 dp dv 2 

\ r dr ) dr + Z dz dz 

Vz dv r / 77 ^ dp\ _ 2 dp .( v T dv r \ 

? + dz \ 3 r + dr ) 3 dz \ r dr ) 


dp . _ d 2 v z 
-Q- z +Jr B ° + T lfo2 + 


Ar, d 2 


r l d 2 Vr 
3 drdz 


d 2 v z 


dv z 


V v t . r.1 . t-. 

d a -2 + q t + r 


dr 2 


! dv ± At, d 2 v z 

! dz 3 dz 2 


= - S z 


( 22 ) 


(23) 
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where 


'rj dr] 
r + d r 


dp . dv T fp dr}\ 2 dr] (v T 0Vr \ q d\r 

dz 9 dz 1 3r dr) 3 dz l r dr ) 3 drdz 


Energy. The electron and ion energy equations are adapted from Mitchner and Kruger, 27 converted to steady- 
state form, and rewritten in terms of the plasma mass density. Grouping terms by order, the electron energy 
equation may be written: 

, BT. , ST, _ , djT, _ 

“• + dr +r - dz r ‘ T ’ + K 'gz‘~ 4 ‘ w 

where T ft is the electron temperature, k c is the Spitzer-Harm electron thermal conductivity: 26 

7h(4'7reo) 2 k D (k B T c y 1/2 m-iu Tl ^ 


Ir- - » 2.55 x 10‘ lo fV 
327r 1 / 2 m 1 ' /2 e 4 In A e lnA e 


and the equation coefficients are: 


1 K e 3 _ 

r - = + 

rj = (“) 

^ _ / v T dv r dv z \ 

T r, = P R (^7 + + Q z ) 

The electron energy source term S e is given by: 

S e = j ■ (E + v x B) - AQ c i (29) 

where the first term on the right represents joule heating of the plasma electrons, and the second term 
represents an energy exchange between the electrons and plasma ions: 

A Q* = 3jjpR(T.-T i )v ei (30) 

with the electron-ion collision frequency given by: 

1 00c w in — 6 P^ n 

v e i — 1.836x10 o/ 2 

mtV 


The steady-state ion energy equation takes the form: 


1 g r 1 + F ' dr 1 dz ‘ 7 ' dz* 


where T* is the ion temperature, k % is the Spitzer-Harm ion thermal conductivity: 


Ki — 1.35 x 10~ 
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and the equation coefficients are: 


r , 1 
r , 2 
r ? 


Ki dKi 3 
_ + _ - - PRVr 

dKi 3 

17 - - 2 pRv - 



dv r dv z \ 
dr^dz) 


The ion energy source term Si is given by: 

Si = $v + AQci 


( 33 ) 


(34) 


where the first term on the right represents ion heating by viscous dissipation, and the second term again 
represents the energy exchange between the electrons and ions. Joule heating of the ions is neglected in the 
ion source term, as it is assumed the ion current density is substantially smaller than the electron current 
density. 

The viscous dissipation $v is : 26 



where A represents the second coefficient of viscosity for a monatomic gas: 

A + -77 - 0 (36) 


The set of fluid equations are closed with an ideal gas equation of state relating the plasma pressure, density, 
and temperature: 

p=( W +Pe) = P«(r<+Te) (37) 


lie. Finite Difference Formulation 


The governing equations have the general form: 


d 2 Y dY dY d 2 Y v - 

f = a ^+ a ^ + a ^ +a ^ +a ° Y + s 


(38) 


where the a n denote nonconstant coefficients and S represents a possible source term. The coupled form 
of the equations hinders an analytic solution, and suggests the use of an iterative numerical approach. The 
equations are written in finite difference form, and solved on a uniform rectangular grid using a generalized 
Newton-Raphson iteration method . 28 Second derivatives are represented by a second-order central difference 
formulation: 


d 2 Y = y(j l i + l)-2y(j,i) + y(j > i- 1) 

dr 2 (Ar ) 2 

d 2 Y = y(j + M)- 2 y(j, i) + y( j-M ) 

dz 2 
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where the index j represents nodal points along the axial grid direction, i represents nodal points along the 
radial grid direction, and A r and A z are the radial and axial grid spacings, respectively. First derivatives 
are written using a first-order switching scheme, which is based upon the sign of the coefficients multiplying 
the derivatives: 28 


dY _ I 

r ULitn-mti 
1 Ar 

a 2 > 0 


dr 1 

I ru,i)-YU,i-D 

( Ar 

c *2 < o 

(40) 

dY I 

, ro+i.o-ro.o 

| A, 

C*3 > 0 


dz ~ | 

I YU,i)-YU-l,i) 

' Az 

013 < 0 



This scheme preserves the dominance of the F(j, z) term and assures stability of the numerical solution. 

Once written in their finite difference analogs, each equation is regrouped into the general form: 

X - + aiY(j, i+ 1) + a 2 Y(j,i- 1) + a 3 Y{j + l,i) + a 4 Y(j - l,i) + S{j,i) (41) 

where the On are nonconstant coefficients, which must be evaluated at each grid location. The generalized 
Newton-Raphson iteration method 

») = Y{j, i)-U (42) 

is then used to iteratively solve each finite difference equation. An over-relaxation factor uj is used to speed 
convergence; its value lies between 1 and 2, and is determined through trial and error for each equation being 
solved. 

Program execution is presented schematically in Figure 2. The fixed grid is divided into 50 radial by 100 
axial nodes, and each equation is iterated first in the axial and then in the radial direction until a relative 
convergence criteria is satisfied (typically to within a tolerance of 10 -4 ). The process is repeated for each 
equation in turn, and the entire loop through the full set of equations is repeated until the exhaust velocity 
and plasma potential have each converged to within 1% of their previous loop values. 

Thrust, Specific Impulse, and Flow Efficiency 

Two complementary methods are used to calculate the total thrust, providing an independent check of 
the thrust value and allowing an estimate to be made of the various force components. The first method 
calculates the total electromagnetic thrust by numerically integrating the axial j x B body force over the 
total current carrying volume, including regions downstream of the thruster exit. Pressure forces generated 
within the thruster are numerically integrated over the inner surfaces of the thrust chamber, and combined 
with the electromagnetic thrust to estimate the total thrust: 

F = f (j x B) z dVtot 4 / pdS tc (43) 

JViot Js ic 

where V to t denotes the total current carrying volume and St c denotes surface regions within the thrust 
chamber. 

A second thrust calculation is performed using; 

F — 77lV ex -j" f (j X B'jzdVoui -f {pcxit Pt«nfc)A c (4^) 

dV aut 
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( 45 ) 


The velocity at the exit plane of the thruster (v ea: »t) is given by: 

_ Et p{i a ' 

Vextt EiPti a ' i ) v *U a ’ i ) 

where the axial index ja denotes the anode exit plane and the radial integration extends from the thruster 
centerline to the inner anode radius. The integration of the electromagnetic body force is performed over the 
current carrying volume external to the thruster ( V ou t )• The pressure force term corresponds to an imbalance 
between the pressure at the anode exit plane (p ex it) an£ l the background gas pressure (ptanfc)j evaluated over 
the thruster exit area A e . Note that the velocity at the thruster exit plane does not in general equal the 
MPD thruster exhaust velocity, due to the additional electromagnetic acceleration which may occur outside 
the thrust chamber and the imbalance of pressure forces at the exit plane. 


Once the total thrust F is calculated, the specific impulse (I 3p ) and plasma flow efficiency (rj/) are found 
via: 



F 2 

71 } “ j2rhP) 


(46) 


where rh is the propellant mass flow rate, g is the acceleration of gravity (9.8 m/s 2 ), and P is the power 
deposited in the plasma, equal to the plasma voltage multiplied by the discharge current. The model does 
not incorporate electrode effects, hence the total discharge voltagel necessary for a prediction of the total 
thruster efficiency cannot be calculated. Because electrode power losses may consume a significant portion 
of the total thruster power, 3 experimentally measured thruster efficiencies will generally be much lower than 
the calculated plasma flow efficiency rjf. Nevertheless, the flow efficiency provides a valuable benchmark 
against which various thruster geometries can be compared. 


lid. Starting Values and Bou ndary C ond itions 

Code input consists of the thruster discharge current, propellant ion mass (amu), propellant mass flow 
rate, geometric boundary conditions, and a host of additional options such as maximum bulk electrode 
temperatures, background gas pressures (to model facility effects), and initial estimates for the plasma 
transport coefficients. Ion temperatures are typically set to 3000 K at electrode and insulator surfaces, and 
electron temperatures are assumed to be continuous into the surfaces ( dT e /dn 0). Vacuum tank pressures 
are set between 1.3xl0~ 2 and 1.3xl0“ 4 Pa (10“ 4 and 1(T 6 torr, respectively). No appreciable difference 
is found in calculated thruster performance using either background pressure value. The inlet velocity uo is 
assumed to be sonic and uniform at the backplate. The inlet density po is assumed to be uniform and is 
calculated from the 1-D continuity equation, m = p 0 Uo^0i where A 0 is the exposed backplate surface area. 
Radial velocities are initially set to zero, and are defined by symmetry to remain zero along the centerline. 
Initial temperatures throughout the calculation, region are set to the bulk electrode temperatures. Initial 
velocities and densities in the calculation region are set to their values at the backplate. A no-slip boundary 
condition on velocity is employed at all insulator and electrode surfaces. 

The electrodes are modeled as equipotential surfaces. Electric fields which enter perpendicular to electrode 
surfaces satisfy the condition dE/dh— 0, where h denotes a unit vector normal to the surface. Electric fields 
lying along electrode surfaces are set to zero, consistent with equipotential surface boundary conditions. 
Insulator surfaces support parallel electric fields, but zero current (due to an assumed infinite dielectric 
resistivity) and zero perpendicular electric fields (due to an assumed infinite relative dielectric permittivity). 
The magnetic stream line ip is set equal to -poJ /(2 tt) along the backplate, zero along the centerline, and is 
assumed to be continuous at both the outer radial grid boundary and and downstream axial grid boundary. 
Setting ^“0 the thruster exit plane would prevent the current from blowing out of the thruster, while 
setting ip-0 at the downstream grid boundary might artificially compress the current blown downstream. 

tThe total discharge voltage is the sum of the plasma and electrode fall voltages. 
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To shorten loop convergence times, the temperature and velocity components are averaged with their pre- 
vious loop values. This damps out minor fluctuations as a solution is approached, allowing more rapid 
loop convergence. Tests performed with and without smoothing yield essentially identical results. Model 
convergence with smoothing is typically obtained within 10-12 loops through the equations, requiring ap- 
proximately 30 CPU minutes on a VAX-9000 computer. The code has not yet been optimized for a CRAY 
computing system. 


III. Results and Discussion 


Comparison with Experiment 


The two-temperature self-field MPD thruster code was tested against the experimental performance of the 
Princeton University extended anode thruster. 29 The extended anode thruster, shown schematically in Figure 
3, consists of a cylindrical 3.2 cm radius copper anode surrounding a 0.95 cm radius thoriated tungsten 
cathode. The anode and cathode lengths are 21.6 cm and 20 cm, respectively. Argon propellant is injected 
through an insulating boron nitride backplate, with either all of the propellant injected near the anode 
radius, or with half of the propellant injected through an annulus near the cathode base and half through 
a ring of 12 evenly spaced ports located 2 cm radially from the thruster centerline. For an argon mass flow 
rate of 6 g/s, the onset of thruster instability was observed to occur at a discharge current of approximately 
21 kA for 50:50 propellant injection, and at approximately 41 kA for all-anode gas injection. 

The present two-temperature code and a previous single-temperature version 25 of the model were used to 
predict the performance of the extended anode MPD thruster. Calculated and experimental thrust values are 
plotted in Figure 4, measured total voltage-current and predicted plasma voltage-current characteristics are 
shown in Figure 5, and predicted flow efficiencies are displayed in Figure 6. Both models predict the thrust 
fairly accurately, but diverge from the experimental values at currents approaching 21 kA (dashed vertical line 
in Figure 4). No steady-state solution was obtained with the two-temperature model for discharge currents 
above this value. The predicted plasma potentials correctly reproduce the trend of the experimentally 
measured voltage-current characteristics, but are consistently lower by approximately 40 volts, reflecting the 
lack of electrode fall calculations in both models. The slightly higher flow efficiencies predicted by the two- 
temperature model result from the more accurate calculation of the transport coefficients, which primarily 
impact the viscous loss terms (through the coefficient of viscosity, Eqn. 16) and the predicted plasma voltage 
drop (via the electrical conductivity, Eqn. 7a). The single temperature model overestimates the viscosity 
coefficient, which depends upon the (typically lower) ion temperature, and underestimates the electrical 
conductivity, which is a function of the (typically higher) electron temperature. For the cases considered, 
the viscous drag is small compared to the total thrust (typically less than 1%), and may be neglected. 
However, the higher electrical conductivity lowers the predicted plasma potential (and consequently the 
power deposited in the plasma), hence for similar thrusts the calculated flow efficiencies are somewhat 
higher. For the two-temperature model, the flow efficiency increases from 40% at 10 kA to around 50% at 21 
kA, at which point steady-state convergence could not be obtained. As discussed earlier, the flow efficiency 
significantly overestimates the total thruster efficiency. The calculated plasma voltage at 20 kA is 17 volts 
(Figure 5), yielding a flow efficiency (Eqn. 46) of nearly 50%; however, the measured discharge voltage is 
approximately 70 volts, yielding a total thruster efficiency of only 15%. Such discrepancies underscore the 
need to develop and incorporate accurate electrode fall models into MPD thruster simulations. 

Both the one-temperature and two-temperature self-field models assume uniform gas injection at the back- 
plate, which corresponds most closely to the 50:50 propellant split used in the experiments. As noted above, 
the onset of thruster instability -was reported at a discharge current of approximately 21 kA using 50:50 
gas injection. The lack of steady-state convergence in the two-temperature model for discharge currents 
exceeding the onset current was not anticipated, and an attempt was made to determine if such an effect 
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occurred when modeling other cylindrical geometry thrusters. Two thruster geometries were selected from 
the extensive set of experimental onset parameters catalogued by Preble. 10 The first thruster had a cathode 
with a radius of 0.95 cm and a length of 5.8 cm, and an anode with an inner radius of 5.1 cm and a length 
of 6.0 cm. The argon propellant was injected with a 50:50 split at a mass flow rate of 6 g/s. The onset of 
thruster instability was experimentally measured at & discharge current of around 19-20 kA, corresponding 
to a J 2 /m value of 6. 0-7.0 xlO 10 A 2 -s/kg. The second thruster consisted of a cathode with a radius of 0.47 
cm and a length of 3.2 cm, and an anode with an inner radius of 2.8 cm and a length of 3.2 cm, operated with 
an argon mass flow rate of 3 g/s. The onset of instability for this thruster occurred at a discharge current 
of approximately 15-16 kA, corresponding to a J 2 /m value of approximately 7.5 - 8.5 x 10 10 A 2 -s/kg. The 
two-temperature model was run for both thruster geometries over a range of J 2 /m values, with the results 
shown in Figure 7. The solid data points denote converged steady-state code solutions, and the open data 
points signify the range of performance variation obtained as the code attempted to find a steady-state 
solution. The dashed vertical lines represent the experimentally determined minimum and maximum values 
of the onset current for the two thrusters, which varied due to the method of propellant injection and slight 
adjustments to the cathode lengths. In general, the model does not converge for onset parameters within 
20% of the experimentally determined values. This is rather remarkable, as the single-fluid equations do 
not support plasma instabilities, and the propellant is already assumed to be fully ionized. Numerically, the 
lack of convergence appears in the set of fluid equations, characterized by oscillations in both the ion and 
electron temperatures and the fluid density values (which are manifested in the calculated plasma potentials 
and exhaust velocities, preventing the numerical convergence criteria from being satisfied). It is possible that 
the plasma conditions leading to thruster onset can be modeled within the single-fluid approximations, but 
where the plasma would dissipate energy via instability processes, the model instead simply oscillates about 
a steady-state solution. The incorporation of iidditional energy dissipation mechanisms, such as higher-order 
ionization or microinstability effects, may thus provide significant insight into the mechanisms responsible 
for the onset of thruster instabilities, 20 ’ 30 ’ 31 and requires further investigation. 


Geometric Scaling A nalysis 

A limited set of computer runs were performed to assess the effect of geometric scale changes on self-field 
MPD thruster performance. Table 1 lists the matrix of cylindrical thruster geometries chosen for comparison. 
The first geometry set (MPDT-1) consisted of a 0.5 cm radius cathode, surrounded by a 2.5 cm radius anode. 
The second set (MPDT-2) doubled the anode radius to 5.0 cm, but kept the cathode radius at 0.5 cm. The 
third set kept the anode radius at 5.0 cm, but increased the cathode radius to 1.0 cm. For each combination 
of radii, the electrode lengths vrere scaled from 1 to 5 times the anode radius (with equal anode and cathode 
lengths assumed), yielding 15 cylindrical thruster geometries. The argon propellant mass flow rate was kept 
constant at 1 g/s, and was assumed to be fully ionized for all cases. Four J 2 /m values, ranging from 2.5 x 10 
A 2 -s/kg to l.OxlO 11 A 2 -s/kg, were evaluated for each geometry. Results of the numerical simulations for 
each thruster geometry are presented below. 


MPD T-1. Thrust characteristics for the 2.5 cm anode radius, 0.5 cm cathode radius thruster are 

presented in Figure 8a. The thrust remains fairly constant with increasing electrode length for J 2 /m values 
of 2.5xl0 lo A 2 -s/kg and 5. Ox 10 lo A 2 -s/kg. At J 2 /m~7.5x 10 lo A 2 -s/kg the thrust remains fairly constant 
for l a / T a (anode length to anode radius) values less than 3; however, steady-state solutions could not be 
obtained for / a /r a > 4. Only one steady-state solution was found for J 2 jrh— l.Ox 10 n A 2 -s/kg, at l a / r n — 1? 
the model would not converge for larger electrode lengths. The maximum thrust of 16.5 N occurs for l a / r a 
and J 2 /m — 1.0 x 10 11 A 2 -s/kg, corresponding to a specific impulse of 1680 seconds. 

If a lack of convergence in the steady-state code does reflect unstable thruster operation, as discussed in 
previous sections, then it would appear for this geometry that higher J 2 /rfi values drive the design toward 
shorter electrode lengths, with optimal l a /r a ratios of around unity. It is thus tempting to conclude that 
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J 0 /r 0 ratios should be kept at unity for this thruster under all operating conditions; however, competing 
with this effect is the requirement for high flow efficiencies. Figure 8b plots the flow efficiency versus l a /r a as 
a function of J 2 jrn, At 2.5 x 10 lo A 2 -s/kg, the flow efficiency dramatically increases with increasing thruster 
length from roughly 48% at / a /r a = 1 to over 60% at l a /r a = 5, where the flow efficiency begins to level 
off. The increasing efficiency is due primarily to a decrease in the calculated plasma potentials (Figure 8c); 
longer electrodes provide a lower current density, which decreases the electric field and lowers the plasma 
potential. At the higher J 2 /rh value of 5.0 x 10 ln A 2 -s/kg, the flow efficiency increases from roughly 63% at 
l a /r a — 1 to around 68% at / a /r a — 2, after which the efficiency remains fairly constant with increasing 
thruster length. At the still higher J 2 jrn value of 7.5 x 10 lo A 2 -s/kg, the flow efficiency only slightly increases 
with increasing electrode length, and the simulation was unable to converge for l a /T a > 4. The efficiency 
at J 2 jrn = 1.0 x 10 11 A 2 -s/kg and / a /r a =1 was approximately 65%, comparable to the efficiency obtained 
at J 2 jrn = 7.5 x 10 lo A 2 -s/kg for the same l a /r a ratio. Steady-state convergence was not possible for 
J 2 jrn — 1.0 x 10 n A 2 -s/kg and l a /r a > 2. As noted previously, a trade-off must thus be made between 
efficient operation at lower J 2 jrn values, which requires longer electrode lengths, and thruster stability 
requirements, which drive the electrode design toward shorter electrode lengths for higher values of J 2 jrn. 

MPDT-2. The second geometry set chosen for simulation retained the same cathode radius of 0.5 cm, 
but doubled the anode radius from 2.5 cm to 5.0 cm for a thruster aspect ratio of 10. The steady-state code 
was able to converge for only a few operating conditions, as illustrated in Figures 9a through 9c. Convergence 
was obtained for J 2 jrn values of 2.5 x 10 lo A 2 -s/kg when l a /r a < 3, and 5. Ox 10 lo A 2 -s/kg when l a /r a = 1. 
Convergence could not be obtained for higher J 2 jrn values, even with the shorter electrode lengths. Thrust 
(Figure 9a) decreases with increasing electrode length at the lower J 2 jrh value, while the flow efficiency 
(Figure 9b) remains constant out to / a /r a = 2, after which the flow efficiency also declines. The maximum 
calculated thrust, at l a /r a — 1 and J 2 /m=5.0x 10 lo A 2 -s/kg, was approximately 16 N, corresponding to 
a specific impulse of roughly 1630 s. These results indicate that, for this larger aspect ratio thruster, 
shorter electrodes are preferable under all operating conditions. The flow efficiency in the stable operating 
regions is approximately 75%, higher than the flow efficiencies found under similar operating conditions in 
MPDT-1. Unfortunately, MPDT-2 appears to be limited to rather low values of J 2 /rh before the onset of 
thruster instability, and because specific impulse scales as J 2 jrn ) the utility of such a device would be greatly 
restricted. The reduced stability of this large aspect ratio thruster agrees with the general results of Preble’s 
onset analysis, 10 which found that smaller aspect ratio thrusters achieved higher values of J 2 jrn before the 
onset of thruster voltage oscillations, 

MPDT-3. The third set of geometries investigated, MPDT-3, retained the larger anode radius of 5.0 
cm but doubled the cathode radius to 1.0 cm, returning the thruster aspect ratio to 5. Predicted thrust, 
voltage, and flow efficiencies versus / a /r a are plotted in Figures 10a through 10c as functions of J 2 jm. The 
reduction in aspect ratio increased the parameter space for steady-state code convergence, indicating that 
stable thruster performance may be obtained over a wider range of operating conditions for small aspect ratio 
thrusters (in agreement with Preble’s onset survey). Thrust values for the lower J 2 jrn values agree fairly well 
with the thrust values predicted for MPDT-1, which is to be expected for self-field thrusters with similar 
aspect ratios. However, only one steady-state solution was obtained for MPDT-3 at J 2 /m=5.0x 10 lo A 2 - 
s/kg, at l a /r a — 1, and no steady-state solution was obtained at J 2 /rh=1.0x 10 11 A 2 -s/kg. Although the 
aspect ratios were the same, the larger electrode radii in MPDT-3 reduced the stable operating regime of 
the thruster. Additional differences are observed in the flow efficiencies of each thruster. Whereas the flow 
efficiency steadily increased with l a /r a at J 2 /m=2.5x 10 lo A 2 -s/kg in MPDT-1, the flow efficiency increased 
and peaked at J a /r a =4 in MPDT-3, after which the flow efficiency decreased with increasing electrode 
length. At J 2 /m=5.0x 10 lo A 2 -s/kg, the flow efficiency increased in MPDT-1 until l a /r a = 2, after which the 
flow efficiency remained approximately constant. However, the flow efficiency steadily decreased in MPDT-3 
with increasing electrode length under the same operating conditions. Both trends may be explained by 
comparing the thrust (Figure 10a) and plasma voltage values (Figure 10c). At J 2 /rn~ 2.5x 10 lo A 2 -s/kg, 
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the thrust remains fairly constant from I 0 /r a = 1 to 4, then drops approximately 5% as l a /r a is increased 
from 4 to 5. The voltage, however, drops almost 30% as l a /r a increases from 1 to 4, then remains fairly 
constant as l a /r a changes from 4 to 5. Referring to Equation 46, the nearly constant thrust, coupled with 
the decreasing plasma voltage, results in an increasing plasma flow efficiency for 1 < l a /r a < 4. As l a /r a 
is increased to 5, the plasma voltage levels off and t he thrust decreases, resulting in a lower plasma flow 
efficiency. At J 2 /m=5.0x 10 lo A 2 -s/kg, both the thrust and plasma voltage are decreasing with increasing 
la/i'a. in such a manner that the plasma flow efficiency decreases with increasing l a /r a as well. It is of 
interest to note that MPDT-1 is the half-scale counterpart of MPDT-3, and in partial agreement with 
Gilland’s results 11 the performance of these half-scale and full-scale thrusters are indeed similar for a limited 
parameter space. However, for these cylindrical self-field thrusters, the similarity in performance starts to 
diverge at higher values of J 2 /m and larger values of l a /r a , even though the mass flow rates were equivalent. 
The maximum calculated thrust for MPDT-3 was roughly 14 N for la/i'a — 1 nnd J /rn— 7.5x10 A -s/kg, 
which corresponds to a maximum specific impulse impulse of approximately 1430 s. 


Stability Diagrams. Based on the preceding results, stability diagrams were generated for the three 
sets of thruster geometries (MPDT-1 through 3) as a function of l a /r a versus J 2 jrn (Figures 11-13). Regions 
of steady-state code convergence are depicted as “stable”, and nonconvergent regimes as “unstable”. Regions 
without data were left unshaded. Though limited in extent, fairly definite patterns emerge from the diagrams. 
In general, the thruster simulations were nonconvergent at higher values of J 2 jrn and larger values of l a /r a , 
with the effects becoming more pronounced for the larger aspect ratio thruster. The smaller diameter 
thruster (MPDT-1) had the largest stable operating space, followed by MPDT-3 (same aspect ratio, but 
twice the size) and MPDT-2 (with the same anode diameter as MPDT-3 but twice the aspect ratio). Of 
particular interest are the slopes of the lines denoting stable operation; the slopes are identical for MPDT-1 
and MPDT-2, which had identical cathode radii bu t dif ferent aspect ratios. The stability region for MPDT-3, 
in which the cathode radius was doubled, has twice the slope of the other stability regions, even though the 
aspect ratios of MPDT-1 and MPDT-3 are identical In addition, the maximum value of J 2 /m for stable 
operation was lower for the larger diameter thruster (MPDT-3) with the same aspect ratio as the smaller 
thruster (MPDT-1), and was halved when the thruster aspect ratio was doubled (MPDT-2 versus MPDT-1). 


These intriguing correlations prompted an attempt to derive a stability scaling relation for the three thruster 
geometries, and a scaling relation was found of the form: 


6.25 x 10 9 / 1 
\rh J c ~ r c \l a 


-fe)-(T)] 


A 2 — s/kg 


(47) 


where (j 2 /? 7 1 ) denotes the maximum J 2 jrn value (in A 2 -s/kg) for stable thruster operation and r tt , l a , r ci 
and l c are the anode radius and length and cathode radius and length, respectively, measured in centimeters. 
The inverse scaling with cathode radius is apparent from the stability diagrams, and the factors of 5 and 
10 appearing in the expression are apparently related to the maximum values of l a j^a an ^ r «/ r c used in 
this limited study. The remaining terms are not intuitively obvious, however, and an effort was undertaken 
to determine the validity of the scaling relation for other cylindrical self-field thrusters. Once again, the 
extensive data base compiled by Preble 10 proved invaluable for this task. A comparison of experimentally 
determined onset values with the values predicted by Eqn. 47 are listed in Table 2 for thruster geometries 
which fall within the geometric constraints of the model: straight, cylindrical self-field thrusters, with uniform 
argon propellant injection, which satisfy the conditions l a /v a < 5, r a jv c < 10, r a > 2.54 cm, r c > 0.5 cm, 
and approximately equal electrode lengths ( l a — l c ). The numerical simulations assumed equal electrode 
lengths, and the factor of (f c // a ) in Eqn. 47 arises from an examination of the Preble data set. In general, 
the predictions of Eqn. 47 and the experimental data agree fairly well, with better conformity at shorter 
length electrodes. Thus, for thruster geometries which fall within the geometric constraints outlined above, 
the stability scaling relationship may be useful as a predictor of thruster instability. 
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The above scaling relation was used to predict the stability performance of an MPD thruster with a cathode 
radiuB of 1 cm and anode radius of 10 cm, the full-scale counterpart to the MPDT-2 thruster above (with 
an aspect ratio of 10). Based on the above results for MPDT-2, a value of l a /r a = 1 was chosen to improve 
the region of operational stability. Inserting values for r a , r c , and / a /r 0 into Eqn. 47 yields a maximum 
J 2 /m value of 2.5 x 10 lo A 2 -s/kg for stable code convergence (and, based on the previous arguments, thruster 
stability). The two-temperature model was then used to simulate the thruster, with a mass flow rate of 1 g/s 
(argon) and discharge currents of 5 kA and 7.5 kA; the model converged to a steady-state solution for the 
5 kA discharge current, but did not converge for the 7.5 kA value. The mass flow was then raised to 4 g/s, 
with discharge currents of 10 kA and 14 kA; the simulation converged for the lower discharge current, but 
no steady-state solution could be obtained for the higher value. The thruster geometry was then changed, 
keeping the anode and cathode radii the same but extending l a /r a to 3. The stability equation predicts 
a maximum J 2 /m value of 1.25 x 10 lo A 2 -s/kg. The code was operated with a mass flow rate of 1 g/s at 
discharge currents of 3.5 kA and 5 kA; the model converged for the lower discharge current but not for the 
higher value. It must be reiterated that the above results only show that the stability equation, based upon 
the performance of a numerical model, can predict the stable regions of model operation; it is the somewhat 
tenuous association between steady-state model convergence and experimentally measured onset data that 
allows such speculation concerning stable thruster operating regimes. 

Based on the above results, and their associated limitations, the following observations may now be tenta- 
tively added to the list of maxims discussed in the introduction. For cylindrical self-field MPD thrusters, 
operated on argon and satisfying the geometric constraints listed above, stable operation at higher values of 
J 2 /m requires smaller ratios of / 0 /r a . At lower values of J 2 /m, longer electrodes are required to increase 
the thruster flow efficiency. Smaller aspect ratio thrusters are generally more stable than larger aspect ratio 
thrusters of equivalent anode diameter, and larger aspect ratio thrusters require shorter electrode lengths 
for stable operation. Smaller diameter thrusters are in general more stable over a wider range of operating 
conditions than their large-scale geometric counterparts. These preliminary results are based on a very lim- 
ited set of numerical data, and such important effects as partial versus full ionization, nonequivalent anode 
and cathode lengths, nonuniform propellant injection, different propellant species, and the effects of flared 
anode geometries (to name but a few) remain to be investigated. 


IV. Concluding Remarks 


A two-dimensional, two-temperature, single fluid magnetohydrodynamics code was developed to predict 
steady-state, self-field MPD thruster performance. The governing equations and numerical methods of 
solution were outlined and discussed. Experimental comparisons with the Princeton extended anode MPD 
thruster were used to benchmark model predictions. The model accurately predicts thrust and reproduces 
trends in the discharge voltage for discharge currents below experimentally measured onset values. However, 
because the model does not include electrode effects the calculated voltage drops are significantly lower 
than experimentally measured values. Predictions of thrust and flow efficiency are made for a matrix of 15 
cylindrical thruster geometries over a range of J 2 /m values, assuming fully ionized argon propellant with 
a mass flow rate of 1 g/s in each case. The simulations indicate that thruster operation at high values of 
J 2 fm requires short electrode lengths for stable operation. At lower values of J 2 /m, longer electrodes are 
required to improve the thruster flow efficiency. Small aspect ratio thrusters are in general more stable than 
larger aspect ratio thrusters of equivalent anode diameter, and larger aspect ratio thrusters require shorter 
electrode lengths for stable operation. Smaller diameter thrusters are generally more stable over a wider 
range of operating conditions than their large-scale geometric counterparts. A maximum specific impulse of 
1680 s was achieved with a 2.5 cm anode radius, 0-5 cm cathode radius thruster with l a /r a = 1 operated 
at J 3 /m = 1.0 x 10 11 A 2 -s/kg. A scaling relation to predict the onset of self-field thruster instability as a 
function of geometry was derived and tested against published onset data. Within the constraints imposed 
by the model, the scaling relation could generally predict onset to within 20% of the experimental values. 
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Table 1. Matrix of thruster geometries and J 2 /m values used in numerical simulations. r OJ r c , and l a denote 
anode radius, cathode radius, and electrode length, respectively. M S” denotes steady-state solutions, “U” 
denotes lack of steady-state code convergence. All data for m = 1 g/s, Ar. 
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Table 2. Comparison of experimentally determined onset values (EXPT) with values predicted by stability 
scaling relation (Equation 47). 









PROGRAM OUTLINE 



Figure 2. Program outline of MPD thruster numerical simulation. 



Figure 3. Schematic of Princeton University extended anode MPD thruster 29 . 
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THRUST CHARACTERISTICS 
EXTENDED ANODE MPDT (Ar: 0 g/s) 


VOLTAGE-CURRENT CHARACTERISTICS 
EXTENDED ANODE MPDT (Ar: 6 g/s) 
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Figure 4. Thrust characteristics for the extended anode MPD 
thruster, m ss 8 g/s (Ar). 1-T: single temperature code, 2-T: 
two- temperature code results. Dashed vertical line corresponds 
to experimentally measured onset current with 50:50 propellant 
injection. 


Figure 5. Measured total voltage (— ) and calculated plasma 

voltage (symbols) for the extended anode MPD thruster, m = 
6 g/s (Ar). 1-T: single temperature code, 2-T: two-temperature 
code results. Dashed vertical line corresponds to experimentally 
measured onset current with 50:50 propellant injection. 


FLOW EFFICIENCY vs. DISCHARGE CURRENT 
EXTENDED ANODE MPDT (Ar: 0 g/s) 



Figure 0. Calculated plasma flow efficiency for the extended an- 
ode MPD thruster, m = 0 g/s (Ar). 1-T: single temperature 

code, 2-T: two-temperature code results. Dashed vertical line 
corresponds to experimentally measured onset current with 50:50 
propellant injection. 


STEADY-STATE CONVERGENCE TEST 
PREBLE THESIS COMPARISON 



Figure 7. Steady-state, two-temperature model convergence tests. 
MPDT-a: r ft = 5.1 cm, r c = 0.95 cm, l a = 6.0 cm, l c = 5.8 cm, m 
= 6 g/s (Ar). MPDT-b: r fl = 2.8 cm, r c = 0.475 cm, = 3.2 cm, 
l c - 3.2 cm, m = 3 g/s (Ar). Dashed vertical lines correspond to 
experimentally measured minimum and maximum onset current 
values. (S): Stable, (U): Unstable code convergence. 
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MPDT-1: THRUST vs. LENGTH/RADIUS 
(r„ = 2.5 cm, t c = 0.5 cm, 1 g/s Ar) 


MPDT-1: FLOW EFFICIENCY v>. I.F.NGTH/KADIUS 
(r„ = 2.5 era, r c = 0.5 cm. 1 g/s Ar) 


J l fm 
O HOP TO 




THRUSTER LENGTH / ANODE RADIUS THRUSTER LENGTH / ANODE RADIUS 


Figure 8(a). Calculated thrust characteristics for MPDT-1. *0>)- Calculated flow efficiencies for MPDT-1 (steady- 

(S): stable, (U): unstable code convergence. state solutions). 


MPDT-1: VOLTAGE vs. LENGTH/RADIUS 
(r. = 2.5 cm, r f = 0.5 cm, I g/s Ar) 



Figure 8(c). Calculated plasma potentials for MPDT-1 (steady- 
state solutions). 
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THRUST (N) 


MPDT-2: THRUST vs, LENGTH/RADIUS 
( r « = 5-0 cin, r c = 0.5 cm, 1 g/s Ar) 


MPDT-2: FLOW EFFICIENCY vs. LENGTH/RADIUS 
(r« = 5.0 era, r c = 0.5 cm, I g/s Ar) 



Figure 9(a)* Calculated thrust characteristics for MPDT-2. Figure 9(b). Calculated flow efficiencies for MPDT-2 (steady- 

(S): stable, (U): unstable code convergence. state solutions). 


MPDT-2: VOLTAGE vs. LENGTH/RADIUS 
(r a = 5.0 cm, r c = 0.5 cm, 1 g/s Ar) 



THRUSTER LENGTH / ANODE RADIUS 


Figure 9(c). Calculated plasma potentials for MPDT-2 (steady- 
state solutions). 


24 


MPDT-3: THRUST vs. LENGTH/RADIUS 
(r. = 5.0 cm, r c = 1.0 cm, I g/s Ar) 


MPDT-3: FLOW EFFICIENCY vs. LENCTII/RADIIS 
(r* = 5.0 cm, r c — 1.0 cm, 1 g/s Ar) 



THRUSTER LENGTH / ANODE RADIUS 



Figure 10(a). Calculated thrust characteristics for MPDT-3 
(S): stable, (U): unstable code convergence. 


Figure 10(b). Calculated flow efficiencies for MPDT-3 (steady 
state solutions). 


MPDT-3: VOLTAGE vs. LENGTH/RADIUS 
(r® = 5,0 cm, r c = 1.0 cm, 1 g/s Ar) 



THRUSTER LENGTH / ANODE RADIUS 


Figure 10(c). Calculated plasma potentials for MPDT-3 (steady- 
state solutions). 
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MPDT-1: STEADY-STATE MODEL CONVERGENCE 



Figure 11. Numerical stability diagram, MPDT-1. Symbols de- 
note numerical data points; “stable” denotes code convergence, 
“unstable” denotes lack of steady-state convergence. Sloped 
lines are interpolations to the data points. 


MPDT-2: STEADY-STATE MODEL CONVERGENCE 



Figure 12. Numerical stability diagram, MPDT-2. Symbols de- 
note numerical data points; “stable” denotes code convergence, 
“unstable” denotes lack of steady-state convergence. Sloped 
lines are interpolations to the data points. 


MPDT-3: STEADY-STATE MODEL CONVERGENCE 



Figure 13. Numerical stability diagram, MPDT-3. Symbols de- 
note numerical data points; “stable” denotes code convergence, 
“unstable” denotes lack of steady-state convergence. Sloped 
lines are interpolations to the data points. 


26 




REPORT DOCUMENTATION PAGE 


Form Approved 
OMB No. 0704 0188 


Public repotting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 
gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of thts 
collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for information Operations and Reports, 1215 Jefferson 
Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503. 


1. AGENCY USE ONLY ( Leave blank) 

2. REPORT DATE 

August 1992 

3. REPORT TYPE AND DATES COVERED 

Final Contractor Report 

4. TITLE AND SUBTITLE 

Numerical Simulation of Geometric Scale Effects in Cylindrical Self-Field 
MPD Thrusters 

5. FUNDING NUMBERS 

WU-506-42-31 

6. AUTHOR (S) 

M. LaPointe 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

Sverdrup Technology, Inc. 

Lewis Research Center Group 
2001 Aerospace Parkway 
Brook Park, Ohio 44142 

8. PERFORMING ORGANIZATION 
REPORT NUMBER 

E-7277 

9. SPONSORING/MONITORING AGENCY NAMES(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135-3191 

10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 

NASA CR- 189224 


11. SUPPLEMENTARY NOTES 

Project Manager, James S. Sovey, NASA Lewis Research Center, (216) 433-7454. Prepared for the 28th Joint 
Propulsion Conference and Exhibit cosponsored by the AIAA, SAE, ASME, and ASEE, Nashville, Tennessee, 
July 6-8, 1992 (work funded by NASA Contract NAS3-25266). 


12a. DISTRIBUTION/AVAILABILITY STATEMENT 


12b. DISTRIBUTION CODE 


Unclassified - Unlimited 
Subject Category 20 


13. ABSTRACT (Maximum 200 words) 

A two-dimensional, two-temperature, single fluid magnetohydrodynamics code which incorporates classical plasma 
transport coefficients and Hall effects has been developed to predict steady-state, self-field MPD thruster performance. 
The governing equations and numerical methods of solution are outlined and discussed. Experimental comparisons are 
used to validate model predictions. The model accurately predicts thrust and reproduces trends in the discharge voltage 
for discharge currents below experimentally measured onset values. However, because the model does not include 
electrode effects the calculated voltage drops are significantly lower than experimentally measured values. Predictions 
of thrust and flow efficiency are made for a matrix of fifteen cylindrical thruster geometries assuming a fully ionized argon 
propellant. A maximum predicted specific impulse of 1680 s is obtained for a thruster with an anode radius of 2.5 cm, a 
cathode radius of 0.5 cm, and equal electrode lengths of 2.5 cm. A scaling relation is developed to predict, within limits, 
the onset of cylindrical, self-field thruster instability as a function of geometry and operating condition. 



17. SECURITY CLASSIFICATION j 18. SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION 20. LIMITATION OF ABSTRACT 
OF REPORT OF THIS PAGE OF ABSTRACT 


Unclassified Unclassified Unclassified 


NSN 7540-01-280-5500 


Standard Form 298 (Rev. 2-89) 
Prescribed by ANSI Std Z39-18 
298-102 













